This notebook shows how to plot and get the nodes of a convex polyhedral function.
Consider the polyhedral function $$ f(x) = \max(-4x - 1, 2x - 1, -x/4, x/2) $$ in the interval $x \in [-1, 1]$.
As the function is convex (as it is the maximum of linear functions), its epigraph is convex. It is defined as the set of $x, y$ satisfying \begin{align*} -1 & \le x\\ x & \le 1\\ -4x - 1 & \le y\\ 2x - 1 & \le y\\ -x/4 & \le y\\ x/2 & \le y \end{align*} We can create this H-representation in Polyhedra as follows:
In [1]:
using Polyhedra
h = HalfSpace([-1, 0], 1) ∩ HalfSpace([1, 0], 1) ∩ HalfSpace([-4, -1], 1) ∩ HalfSpace([2, -1], 1) ∩
HalfSpace([-1/4, -1], 0) ∩ HalfSpace([1/2, -1], 0)
Out[1]:
We now transform it to a polyhedron using CDDLib, the library that is chosen to compute the V-representation.
In [2]:
using CDDLib
p = polyhedron(h, CDDLib.Library())
Out[2]:
Note that creating the polyhedron does not trigger the computation of the V-representation by itself.
This is done when the V-representation is queried, e.g. by vrep.
We can see that 5 nodes of the polyhedral function with the ray $(0, 1)$ which is expected for an epigraph.
In [3]:
vrep(p)
Out[3]:
We see now that the V-representation is known by the polyhedron. Now both the H- and V-representation are available.
In [4]:
p
Out[4]:
In [5]:
using Plots
plot(p ∩ HalfSpace([0, 1], 4))
scatter!([x[1] for x in points(p)], [x[2] for x in points(p)])
Out[5]:
JuMP variables can be constrained to be in the epigraph.
Note that if the optimization model simply involves minimizing or maximizing a linear objective function over the polyhedron, the problem can simply the solved by iterating over the V-representation.
One can use Polyhedra.linear_objective_solver to get VRepOptimizer when the V-representation is computed.
This optimizer does exactly what we just described.
In [6]:
using JuMP
model = Model(Polyhedra.linear_objective_solver(p))
@variable(model, x)
@variable(model, y)
@constraint(model, [x, y] in p)
@objective(model, Min, x + y)
optimize!(model)
The optimal solution of this problem is one of the extreme points of p:
In [7]:
termination_status(model)
Out[7]:
In [8]:
value(x), value(y)
Out[8]:
In [9]:
plot!([-0.2, -1.2], [0.0, 1.0], linewidth=2) # Objective function
scatter!([value(x)], [value(y)], markershape=:star5, markersize=8) # Optimal solution
Out[9]:
In [10]:
@objective(model, Max, x + y)
optimize!(model)
As we can see, maximizing $x + y$ is an unbounded problem, the infeasibility certificate is given by the extreme ray $(0, 1)$.
In [11]:
termination_status(model)
Out[11]:
In [12]:
value(x), value(y)
Out[12]:
The epigraph can also be used in a more complex optimization problem involving other variables and constraints. Here Polyhedra.linear_objective_solver should not be used because the feasible set of the new model is not exactly p so the V-representation will have to be recomputed which is less efficient than solving the problem using a linear programming solver such as GLPK.
In [13]:
using JuMP
import GLPK
model = Model(with_optimizer(GLPK.Optimizer))
@variable(model, x <= -2)
@variable(model, y)
@variable(model, 0 <= u <= 3)
@constraint(model, [2x + u, y] in p)
@objective(model, Min, u + y)
optimize!(model)
In [14]:
termination_status(model)
Out[14]:
In [15]:
value(x), value(y), value(u)
Out[15]:
Consider the polyhedral function $$ f(x) = \max(-4x - 2y - 1, 2x - y - 1, -x/4 + y/2, x/2 + y) $$ in the interval $(x, y) \in [-1, 1]^2$.
As the function is convex (as it is the maximum of linear functions), its epigraph is convex. It is defined as the set of $x, y, z$ satisfying \begin{align*} -1 & \le x\\ x & \le 1\\ -1 & \le y\\ y & \le 1\\ -4x - 2y - 1 & \le z\\ 2x - y - 1 & \le z\\ -x/4 + y/2 & \le z\\ x/2 + y & \le z \end{align*} We can create this H-representation in Polyhedra as follows:
In [16]:
using Polyhedra
h = HalfSpace([-1, 0, 0], 1) ∩ HalfSpace([1, 0, 0], 1) ∩
HalfSpace([0, -1, 0], 1) ∩ HalfSpace([0, 1, 0], 1) ∩
HalfSpace([-4, -2, -1], 1) ∩ HalfSpace([2, -1, -1], 1) ∩
HalfSpace([-1/4, 1/2, -1], 0) ∩ HalfSpace([1/2, 1, -1], 0)
Out[16]:
We create the polyhedron the same way as for 1-dimensional polyhedral functions.
In [17]:
using CDDLib
p = polyhedron(h, CDDLib.Library())
Out[17]:
We obtain the 10 nodes as follows along with the $(0, 0, 1)$ ray that is expected as the polyhedron is an epigraph.
In [18]:
vrep(p)
Out[18]:
In [ ]:
m = Polyhedra.Mesh(p)
In [ ]:
using MeshCat
vis = Visualizer()
setobject!(vis, m)
IJuliaCell(vis)
In [ ]:
using Makie
Makie.mesh(m, color=:blue)
In [ ]:
using Makie
Makie.wireframe(m)
In [ ]: